### RPS PUBLIC OPINION
### May 16, 2017


library(foreign)
library(dplyr)
library(ggplot2)
library(grid)
library(gtable)
library(foreign)
library(reshape2)
library(stringr)
library(GGally)
library(Hmisc)
library(survey)
library(stargazer)

rm(list=ls()) #to remove all files 

options(digits=10)
options(scipen=10)


load(file="cces_energy_replication.RData")
load(file="state_sample_size.RData")
states<-read.csv(file="StateCodes.csv")

data$cces2014_staterps<-as.numeric(as.vector(data$cces2014_staterps))
library(arm)
multilevel.reg<-(lmer(cces2014_staterps~(1|abb)+scale(log(wind_potential+.01))+scale(log(solar_potential+.01))+factor(education)+factor(pid3)+scale(log(oil_coal_percent+.01))+scale(log(mining_percent+.01)), data=data))
stargazer(multilevel.reg)
survey.data <- svydesign(ids = data$respondent, data = data, weights = data$weight)

## raw mean
prop.table(table(data$cces2014_staterps))
## weighted mean
prop.table(svytable(~data$cces2014_staterps, design = survey.data))

data.graph.weighted<-prop.table(prop.table(svytable(~data$state+data$cces2014_staterps, design = survey.data)),1)[,2]

data.graph<-prop.table(table(data$state, data$cces2014_staterps),1)[,2]
data.graph<-cbind(names(data.graph), as.vector(data.graph))
data.graph<-cbind(data.graph,data.graph.weighted )
colnames(data.graph)<-c("state", "cces2014_rps_raw", "cces2014_rps_weighted")
data.graph<-as.data.frame(data.graph)
data.graph$cces2014_rps_raw<-as.numeric(as.vector(data.graph$cces2014_rps_raw))
data.graph$cces2014_rps_weighted<-as.numeric(as.vector(data.graph$cces2014_rps_weighted))
cor(data.graph$cces2014_rps_weighted, data.graph$cces2014_rps_raw, use="complete.obs")
data.graph<-merge(data.graph, states, by.x="state", by.y="Name", all.x=T, all.y=T)
data.graph<-merge(data.graph, sample_size, by.x="POAbrv", by.y="abb")

rps_policies<-read.csv(file="states_RPS_categorized.csv")
data.graph<-merge(data.graph, rps_policies, by.x="POAbrv",by.y="state", all.x=T)
dim(data.graph)

data.graph<-rbind(NA,data.graph)

data.graph$policy_25<-as.vector(data.graph$policy_25)
dim(data.graph)
data.graph$policy_25[1]<-"Current RPS Policy"

#data.graph<-subset(data.graph, !is.na(data.graph$cces2014_rps_weighted))
data.graph$congruent<-0
data.graph$congruent[(data.graph$policy_25=="binding less than 25%" | data.graph$policy_25=="binding 25% or more") & data.graph$cces2014_rps_weighted>.5]<-1
data.graph$congruent[(data.graph$policy_25=="voluntary" | data.graph$policy_25=="no target") & data.graph$cces2014_rps_weighted<.5]<-1
table(data.graph$congruent)
data.graph<-data.graph[order(data.graph$cces2014_rps_weighted),]

data.order<-order(as.numeric(data.graph$cces2014_rps_weighted), decreasing = FALSE)
data.graph<-data.graph[data.order,]
right.label<-as.vector(data.graph$policy_25)
left.label<-as.vector(data.graph$state)

pdf(paste("rps_states_cces.pdf",sep=""),width=7.5, height=8.5)
 par(mfrow=c(1,1),mar=c(3,6,3,8))
#right.label<-paste("(", round(data.graph$cces2014_rps_weighted*100,0), "%)", sep="")
y.axis<-c(1: dim(data.graph)[1])
plot(data.graph$cces2014_rps_weighted*100, y.axis, axes=F, xlim=c(0,100),ylim=c(.1,dim(data.graph)[1]),pch=19,main="",xlab='',ylab='',type="n")
for (i in 1: dim(data.graph)[1]-1){
	 segments(0, i, 100, i, col="grey")
}
segments(0,0,0,(dim(data.graph)[1]), col="light grey")
segments(50, 0, 50, (dim(data.graph)[1]), col="light grey")

segments(data.graph$cces2014_rps_weighted*100-data.graph$se*100, y.axis, data.graph$cces2014_rps_weighted*100+data.graph$se*100, y.axis, col="black",lwd=2)
points(data.graph$cces2014_rps_weighted*100, y.axis, col="black", pch=19, cex=.75)
axis(1,  pos=0,cex.axis=.8,at = seq(0,100,10))
axis(2, at = c(0,y.axis, (dim(data.graph)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=0,cex.axis=.75, tick=FALSE)
axis(4,at = c(0,y.axis, (dim(data.graph)[1]+2)), pos=100,label = c('',right.label,''), tick=FALSE,las = 2,cex.axis=.75, font=1)
title(main="Support for RPS by State",xlab = "", ylab = "", cex.lab=1)
title(xlab = "Percentage of Public that Supports RPS",cex.lab=1, line=1)
dev.off()

